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1.  Introduction 


Neuroscientists  have  shown  increasing  interest  in  recent  years  in  examining  the  integrated 
activity  of  the  brain’s  circuitry  across  large-scale  networks  (Spoms  et  al.  2004;  Friston  2011; 
Alivisatos  et  al.  2012).  This  idea  is  based  on  the  notion  that  functional  connections  are  formed 
through  synchronously  firing  assemblies  (populations)  of  neurons  (e.g.,  Engel  et  al.  2001;  Varela 
et  al.  2001),  and  experimentalists  have  identified  several  frequency  bands,  such  as  alpha  and 
gamma,  that  relate  changes  in  ongoing  brain  oscillations  to  behavioral  performance  (Klimesch 
et  al.  2007;  Jensen  et  al.  2012).  These  functional  connections  are  dynamically  changing,  and 
within  short  time  periods,  connections  between  neural  populations  can  be  created,  broken,  and 
reestablished  with  other  populations  in  response  to  different  sensory  stimuli  and  different 
processing  demands  (Battaglia  et  al.  2012;  Buschman  et  al.  2012).  Given  the  importance  of  these 
functional  networks  for  communication  across  brain  regions,  electroencephalogram  (EEG) 
experiments  have  increasingly  focused  on  linking  variability  in  functional  networks  to 
performance  on  behavioral  tasks  (Supp  et  al.  2005;  Philiastides  and  Sajda  2006;  Sun  et  al.  2009; 
Sakkalis  2011;  loannides  et  al.  2012;  Blinowska  et  al.  2013). 

Researchers  have  attempted  to  devise  different  approaches  and  methods  for  measuring  functional 
networks  accurately  in  EEG  time  series  that  are  robust  to  noise.  One  class  of  functional  measures 
known  as  effective  connectivity  refers  to  directed  interactions  between  different  brain  regions  or 
cell  populations,  and  it  is  closely  related  to  the  concept  of  generalized  synchronization  (Schiff 
et  al.  1996;  Battaglia  et  al.  2012).  There  are  a  variety  of  signal  analysis  measures  used  for  trying 
to  detect  effective  connections  between  different  neural  populations,  including  Kullback-Eeibler 
Divergence,  Transfer  Entropy,  Partial  Directed  Coherence,  Directed  Transfer  Eunction  (DTE), 
and  the  direct  DTE  (for  review,  see  Blinowska-Cieslak  and  Zygierewicz  2012;  Seghouane  and 
Amari  2012).  In  this  analysis,  we  focus  on  DTE,  which  is  related  to  the  concept  of  Granger 
Causality  and  derived  from  multivariate  autoregressive  modeling  (MVAR)  of  the  time  series. 
Granger- Geweke  causality  was  originally  introduced  in  economic  time  series  analysis  and  used 
to  study  dependencies  in  bivariate  signals  (Granger  1969;  Geweke  1982).  If  incorporating  the 
components  of  the  model  of  the  first  signal  helps  reduce  the  error  of  prediction  in  the  second 
signal,  the  first  signal  is  said  to  have  caused  the  second  signal.  Thus  these  effective  connectivity 
measures  hold  promise  for  capturing  the  temporal  dynamics  of  functional  neural  networks 
(Blinowska  et  al.  2013). 

This  study  investigates  how  well  DTE,  an  effective  connectivity  measure,  captures  the 
underlying  network  connectivity  of  a  simulated,  EEG-like  signal  generated  from  several 
thousand  nonlinear  conductance-based  neuron  models.  We  simulate  6  different  networks,  each 
with  4  interacting  nodes  that  oscillate  within  2  frequency  bands  of  interest  in  experimental  EEG 
studies,  alpha  and  gamma  (Nunez  and  Srinivasan  2006).  Eour  of  the  networks  are  simple 
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feed-forward  networks  with  one  node  driving  3  other  nodes,  and  we  manipulate  whether  the 
driving  or  receiving  nodes  oscillate  at  alpha  or  gamma  frequencies.  Two  of  the  networks  have  a 
more  complex  connectivity  pattern  and  include  feedback  connections  and  a  ring  structure. 

Our  evaluation  of  these  measures  differs  from  the  previous  studies,  which  have  used 
multichannel  signals  generated  by  MVAR  models  with  defined  dependencies  in  the  coefficients 
(Baccala  and  Sameshima  2001)  or  periodic  signals  with  introduced  delays  between  channels  as  a 
means  of  validating  the  measures  (Blinowska  et  al.  2013).  In  other  studies,  bivariate  connectivity 
measures  were  validated  using  neural  mass  models  (David  and  Friston  2003;  Wendling  et  al. 
2009).  In  the  present  study,  we  investigate  the  multichannel  directional  connectivity  measure 
related  to  Granger  causality  applied  to  signals  generated  from  4-node  networks.  Since  the 
simulated  signals  are  created  as  a  superposition  of  post-synaptic  potentials  of  conductance-based 
models  of  neurons,  our  model  better  resembles  an  experimental  EEG  signal,  which  is  created  by 
the  activity  of  firing  neurons. 

Erom  the  results,  we  propose  2  guiding  principles  for  the  application  of  effective  connectivity 
measures  to  experimental  data  to  minimize  spurious  connections:  computing  the  measure  only 
when  there  is  overlapping  power  in  the  spectrum  and  applying  a  conservative  statistical 
threshold.  Interpretation  of  these  connectivity  results  can  be  further  enhanced  by  incorporating 
known  anatomical  constraints. 


2.  Methods 


2.1  Neural  Networks  Generating  EEG-like  Signals 

The  simulated  neural  signals  were  generated  using  a  modified  version  (Suffczynski  et  al.  2014) 
of  the  cortical  network  model  originally  developed  to  study  epileptic  phenomena  (Kudela  et  al. 
1997;  Anderson  et  al.  2012).  The  modifications  included  changes  to  simulate  normal  cortical 
activity  without  inducing  seizures.  This  implementation  integrates  nonlinear  conductance-based 
neuron  model  equations  using  Euler's  method  with  a  fixed  step  size  of  0.01  ms. 

The  model  networks  we  examined  all  consist  of  4  interacting  nodes,  each  node  consisting  of  a 
2-layer  model  of  400  excitatory  neurons  (pyramidal  cells)  and  100  inhibitory  interneurons.  A 
schematic  of  the  arrangement  as  well  as  connectivity  between  layers  is  illustrated  in  Eig.  1 .  The 
current  balance  equation  for  the  excitatory  neurons  has  the  form 

^  ~hyn  ~  ha  ~  ha  ~  h  ~  —  h  ~  4-  (1) 

The  form  of  the  equations  for  the  inhibitory  interneurons  is  identical  to  the  excitatory  neurons 
except  for  the  deletion  of  the  calcium-dependent  potassium  current  (Suffczynski  et  al.  2014).  The 
detailed  equations  and  symbols  definitions  are  provided  in  Appendix  A. 
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Fig.  1  A  simulated  node  consists  of  400 

excitatory  neurons  and  100  inhibitory 
neurons.  Connectivity  between  and  within 
layers  were  randomly  chosen.  Example 
consists  of  one  neuron  in  one  layer 
connected  to  3  neurons  in  another  layer. 

For  each  node,  connections  among  the  500  neurons  were  chosen  randomly  from  a  uniform 
distribution  that  ensures  connections  between  all  neurons  were  equally  likely.  Each  connection 
was  one  of  4  possible  types:  excitatory  to  excitatory,  excitatory  to  inhibitory,  inhibitory  to 
inhibitory,  and  inhibitory  to  excitatory.  The  connection  strength  and  delay  between  the  neurons 
was  also  randomly  selected  from  a  Gaussian  distribution.  In  addition  to  these  connections  within 
each  node,  a  total  of  120  connections  were  randomly  selected  between  each  node  pair  to  create 
the  small  networks  used  in  the  simulations.  The  model  equations  and  connectivity  parameters  are 
provided  in  Appendix  A. 

To  generate  different  frequency  oscillations  in  the  simulated  networks,  we  set  different  levels  of 
stimulation  to  different  nodes.  The  stimulation  to  each  node  was  generated  by  a  Poisson  process 
of  action  potentials  to  the  excitatory  cells.  By  varying  the  rate  parameter  X  of  the  Poisson 
process,  we  can  shift  the  frequency  of  the  node's  synchronous  oscillation.  Model  parameters 
were  chosen  so  that  their  power  spectra  peaked  within  2  experimentally  relevant  frequency  bands 
(Nunez  and  Srinivasan  2006):  approximately  1 1  Hz  to  represent  the  alpha  band  and 
approximately  40  Hz  to  represent  the  gamma  band. 
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The  resulting  oscillatory  signal  from  each  node  was  modeled  as  an  EEG-like  signal  by  summing 
and  low-pass  filtering  the  postsynaptic  potentials  of  all  excitatory  synapses  (Nunez  and 
Srinivasan  2006;  Suffczynski  et  al.  2014).  Eor  each  node,  5  s  of  simulated  signals  are  generated 
from  summation  of  the  postsynaptic  potentials,  and  then  they  are  filtered  using  a  low-pass 
Butterworth  filter  of  order  10  with  cut-off  at  100  Hz  and  decimated.  The  first  and  last  50  ms  are 
discarded  in  the  filtering  process.  This  procedure  resulted  in  simulated  signal  with  a  sampling 
frequency  of  200  Hz.  The  resulting  signal  used  in  the  effective  connectivity  analysis  consisted  of 
980  points.  This  length  of  the  signal  window  also  fulfills  constraints  for  computing  an  MVAR 
model,  ensuring  that  the  number  of  estimated  coefficients  is  less  than  1/lOth  of  the  total  number 
of  points  (Blinowska-Cieslak  and  Zygierewicz  2012). 

We  used  this  procedure  to  generate  6  simulated  networks  that  were  either  simple  feed-forward 
networks  or  more  complex  feedback  networks.  Each  network  had  4  nodes.  The  4  simple 
networks  have  one  node  projecting  to  3  other  nodes,  and  we  vary  whether  the  sending  or 
receiving  nodes  oscillate  within  alpha  or  gamma  frequency  bands.  The  2  complex  networks  have 
feedback  connections  and  a  ring  structure  among  the  nodes,  and  we  again  vary  whether  the 
sending  and  receiving  nodes  oscillate  at  alpha  or  gamma. 

2.2  Effective  Connectivity  Measure 

The  time  series  from  6  simulated  EEG-like  networks  were  analyzed  separately  using  an  effective 
connectivity  measure,  DTE.  Related  to  Granger- Geweke  causality  (Granger  1969;  Geweke 
1982),  this  measure  is  based  on  MVAR  modeling  of  time  series  and  have  been  interpreted  as  a 
multivariate  extension  of  the  Granger  concept  (Baccala  and  Sameshima  2001;  Korzeniewska 
et  al.  2003;  Blinowska  2011;  Blinowska-Cieslak  and  Zygierewicz  2012;  Haufe  et  al.  2013). 

When  the  time  series  data  are  fit  to  an  autoregressive  model  by  minimizing  the  error  terms,  the 
model  coefficients  provide  information  about  time  lag  influences  between  the  signals  and  thus 
capture  the  connections  producing  those  signals.  We  examine  how  well  DTE  captures  the 
underlying  connectivity  among  the  nodes  in  our  simulated  EEG-like  networks. 

The  simulated  signals  for  each  network  were  separately  modeled  as  an  MVAR  process;  that  is, 
we  assume  that  the  signals  can  be  described  as 

x(n)  =  —  k)  +  w(n).  (2) 

In  this  equation,  p  is  the  model  order  and  w(n)  is  interpreted  as  a  white  noise  process  or  the  error 
in  the  prediction.  We  solve  for  the  matrix  model  coefficients  A(k)  by  noting  that  x{n  -  k)  is 
uncorrelated  to  w(n).  This  allows  us  to  write  the  equations  as  a  block-Toeplitz  matrix  problem, 
which  can  be  solved  by  the  Eevinson-Wiggins-Robinson  (EWR)  algorithm.  More  details  on  this 
derivation  or  the  EWR  algorithm  can  be  found  in  Golub  and  Eoan  (2012). 

Connectivity  is  then  calculated  from  the  Eourier  transform  of  the  model  coefficients.  DTE  values 
are  obtained  from  the  transfer  matrix  H(f),  which  is  obtained  by  taking  the  inverse  of  the  Eourier 
transform  of  the  model  coefficients  A(f): 
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dtp..  =  l^0-(/)l  ^  (3) 

where 

Hijif)  ^  A\j(in-\  (4) 

and 

5(/)  =  H(/)VH*(f),  (5) 

where  V  is  the  residual  covariance  matrix  of  noise  process  w(n),  S(f)  is  the  spectral  matrix  of  the 
signal,  and  H*(/)  is  a  Hermitian  transpose  of  matrix  H(/). 

To  accurately  compute  an  MVAR  model,  we  must  determine  the  optimal  model  order.  That  is, 
we  want  to  know  how  many  time-lag  terms  to  include  in  the  MVAR  model.  Too  few  terms  will 
result  in  a  poor  characterization  of  the  original  time  series,  and  the  associated  power  spectrum 
profile  will  be  smooth  and  featureless.  Too  many  terms  results  in  over- fitting  and  can  lead  to 
spurious  peaks  in  the  power  spectrum.  The  model  order  used  in  our  analyses  was  determined  by 
computing  information-theoretic  criteria  and  by  comparing  the  parametric-computed  power 
spectrum  with  the  nonparametric-computed  spectrum. 

Several  information-theoretic  criteria  were  computed  since  they  are  based  on  different 
assumptions  about  the  signal  and  use  different  methods  to  penalize  the  model  for  additional 
parameters  (see  Appendix  B),  including  the  Akaike  Information  Criterion  (AIC),  Hannan-Quin 
Criterion  (HQIC),  Bayesian  Information  Criterion  (BIC),  and  Forward  Prediction  Error  (FPE) 
(Akaike  1974;  Akaike  1981;  Stoica  and  Selen  2004;  Blinowska-Cieslak  and  Zygierewicz  2012). 
Additionally,  the  versions  of  the  AIC,  BIC,  and  HQIC  corrected  for  small  sample  sizes  were 
computed:  AICc,  BICc,  and  HGIc.  These  7  information-theoretic  criteria  were  computed 
separately  for  each  of  the  6  network  topologies  examined  in  this  study,  and  the  results  are  shown 
in  Fig.  2.  The  model  order  is  determined  by  searching  for  a  local  minimum  in  the  plots  across  the 
criteria.  Due  to  different  assumptions  about  the  signal,  some  of  them  may  not  show  the 
minimum,  so  the  model  order  is  selected  from  those  plots  that  do  have  a  minimum.  Here,  the 
minimum  of  the  BIC,  FPE,  and  adjusted  BICc  is  approximately  the  same  for  all  network 
topologies,  so  we  chose  a  model  order  ofp=  10.  Each  time  lag  is  one  sample  of  the  time  series 
data,  so  the  MVAR  model  for  each  simulation  examined  whether  or  not  the  current  signal  can  be 
predicted  from  the  previous  10  samples  from  all  other  nodes.  The  appropriateness  of  p  =  10  for 
these  data  was  also  examined  by  comparing  the  power  spectrum  computed  using  the  auto¬ 
regressive  model  with  one  using  Welch's  method,  and  we  obtained  good  agreement  between  the 
parametric  5(/)and  nonparametric  (Matlab's  pwelch  function)  spectrum  (Brovelli  et  al.  2004). 
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Fig.  2  Values  of  7  information-theoretic  criteria  plotted  against  the  model  order  used.  The  legend  is  labeled 
according  to  the  figures  in  which  the  power  spectra  and  associated  DTF  results  for  the  simulated 
networks  are  presented. 

In  addition  to  computing  the  connectivity  measure  across  all  frequencies,  we  also  integrated  DTF 
only  in  the  frequency  ranges  where  the  power  spectrum  shows  peaks  in  our  2  bands  of  interest, 
alpha  (10-12  Hz)  and  gamma  (35-45  Hz).  Theses  ranges  were  determined  by  the  full  width  at 
half-maximum  value  for  the  peak  in  power  spectrum  in  the  node  where  the  peak  in  an  alpha  or 
gamma  range  was  the  highest. 

2.3  Computing  Statistical  Significance 

We  determined  statistical  significance  of  the  computed  measures  by  shuffling  the  phase  of  the 
signals  using  Theiler’s  method  (Prichard  and  Theiler  1994)  to  generate  surrogate  data.  The 
advantage  of  the  Theiler  method  is  that  it  can  be  easily  applied  to  experimentally  obtained  time 
series,  our  targeted  use  of  this  method  in  future  research.  The  method  takes  the  Fourier  transform 
X(f)  of  the  time  series  x(t)  and  multiplies  each  point  by  a  random  phase.  We  symmetrize  X(f) 
and  take  the  inverse  Fourier  transform,  which  results  in  signals  that  have  the  same  power 
spectrum  as  the  original  signal  but  with  no  correlation  in  the  phases  (Blinowska-Cieslak  and 
Zygierewicz  2012).  This  process  was  repeated  500  times  on  the  signal  obtained  from  the 
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connected  network  simulation,  and  we  generated  a  histogram  of  the  results.  From  this,  we 
obtained  the  distribution  of  the  DTF  values  for  each  frequency.  We  computed  the  95th  percentile 
for  the  distribution  at  each  frequency  point  resulting  in  a  significance  level  of  p  <  0.05. 


3.  Results 


We  computed  connectivity  values  for  each  of  our  6  simulated  networks  using  DTF.  In  each 
figure,  section  A  shows  the  network  topology,  section  B  shows  the  power  spectrum  Sn(f) 
between  0  and  100  Hz  as  computed  from  the  MVAR  model,  and  section  C  shows  the  computed 
DTF.  The  rows  in  section  C  correspond  to  the  sending  nodes  whereas  the  columns  correspond  to 
the  receiving  nodes.  In  each  DTF  plot,  the  black  curves  correspond  to  the  significance  level 
computed  using  Theiler's  method  as  described  in  Section  2.3  of  this  report,  and  the  red  curves 
above  the  black  lines  correspond  to  frequency  values  with  DTF  values  that  were  statistically 
significant  based  on  Theiler’s  method. 

3.1  Simple  Networks 

All  of  these  networks  consisted  of  a  central  driving  node  that  projects  to  3  other  receiving  nodes. 
The  parameters  of  the  nodes  were  tuned  so  that  the  nodes  would  intrinsically  oscillate  either  in 
the  alpha  band  centered  around  13  Hz  or  in  the  gamma  range  centered  around  40  Hz.  We 
examined  the  following  4  cases: 

•  An  alpha  frequency  node  drives  3  alpha  frequency  nodes. 

•  A  gamma  frequency  node  drives  3  gamma  frequency  nodes. 

•  An  alpha  frequency  node  drives  3  gamma  frequency  nodes. 

•  A  gamma  frequency  node  drives  3  alpha  frequency  nodes. 

In  the  first  simulation,  an  alpha  frequency  node  drives  3  alpha  frequency  nodes,  as  shown  in 
Fig.  3A,  where  node  1  acts  as  the  driver  of  the  other  3  nodes  with  simulated  synaptic  connections 
projecting  to  nodes  2-4.  The  power  spectrum  in  Fig.  3B  shows  that  node  1  has  the  most  power 
but  all  nodes  peak  near  13  Hz.  Figure  3C  is  the  computed  DTF  between  0  and  100  Hz.  In  the  top 
row  of  Fig.  3C  are  the  DTF  components  corresponding  to  connections  I  A  2,  I  A  3,  and  I  A  4. 
Given  the  underlying  topology  of  the  network,  one  would  expect  that  only  connections  1  A  2,  3, 
4  would  result  in  significant  values  of  the  DTF.  However,  we  see  that  in  8  of  the  9  other  panels 
in  this  plot,  the  red  DTF  values  are  above  the  significance  threshold  shown  in  black.  In 
particular,  we  see  that  the  DTF  is  significant  in  components  4  A  1  and  4  ->  3  across  all 
frequencies  except  for  a  narrow  band  corresponding  to  the  peak  of  the  power  spectrum  (peaked 
near  14  Hz). 
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Fig.  3  Alpha  node  1  driving  3  alpha  nodes.  A:  connectivity  structure  of  a  simulated 
network.  B:  power  spectrum  (0-100  Hz)  computed  from  the  MVAR  model.  C: 
computed  DTF  values  plotted  against  frequency  (0-100  Hz),  where  the  rows 
correspond  to  the  sending  nodes  and  the  columns  correspond  to  the  receiving  nodes. 
In  each  DTF  plot,  the  red  curves  above  the  black  lines  correspond  to  frequency 
values  with  DTF  values  statistically  significant. 


The  second  simulation  shown  in  Fig.  4  is  identical  to  Fig.  3  except  all  nodes  are  now  generating 
a  signal  in  the  gamma  range.  Similar  to  the  first  simulation,  node  1  has  much  more  power  than 
the  other  3  nodes,  indicating  a  larger  degree  of  synchrony.  Examining  the  top  row  of  Fig.  4C,  we 
see  that  the  components  of  the  DTF  corresponding  to  connections  1  ^  2,  3,  4  are  substantially 
above  the  level  of  significance  for  a  wide  range  of  frequencies,  much  more  so  than  the  same 
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components  seen  in  Fig.  3C.  This  being  said,  there  are  still  spurious  results  that  indieate 
connections  3^1  and  2  ^  3.  If  we  examine  the  component  of  the  DTF  corresponding  to  the 
conneetion  1  ^  3,  we  see  that  the  measure  is  above  signifieanee  over  a  broad  range  of 
frequeneies.  In  fact,  it  is  significant  over  a  range  of  frequeneies  in  which  there  is  no  power  in  the 
corresponding  spectrum.  It  is  impossible  for  there  to  be  an  influence  from  one  region  to  another 
within  this  frequency  band,  and  thus  these  results  show  that  DTF  deteets  spurious  eonneetions. 
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Fig.  4  Gamma  node  1  driving  3  gamma  nodes.  A:  connectivity  structure  of  a  simulated  network. 
B;  power  spectrum  (0-100  Hz)  computed  from  the  MVAR  model.  C:  computed  DTF 
values  plotted  against  frequency  (0-100  Hz),  where  the  rows  correspond  to  the  sending 
nodes  and  the  columns  correspond  to  the  receiving  nodes.  In  each  DTF  plot,  the  red 
curves  above  the  black  lines  correspond  to  frequency  values  with  DTF  values  statistically 
significant. 
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Next  we  examined  a  situation  where  the  driving  node  had  little  effect  on  the  receiving  nodes. 
Figure  5  is  a  simple  network  in  which  the  driving  node  oscillates  in  the  alpha  frequency  band  and 
the  receiving  nodes  oscillate  in  the  gamma  frequency  band.  Importantly,  in  the  power  spectrum 
plot  shown  in  Fig.  5B,  there  is  little  overlap  between  the  spectra  of  the  sending  and  receiving 
nodes.  Examining  the  1  ^  2,  3,  4  components  of  the  DTF  in  Fig.  5C,  the  DTF  values  are  not 
above  significance  in  the  frequency  band  of  the  driving  node.  The  DTF  is  slightly  above 
significance  in  the  1^2  and  1^3  components  near  50  Hz. 


A.  B. 


Fig.  5  Alpha  node  1  driving  3  gamma  nodes.  A:  connectivity  structure  of  a  simulated 

network.  B;  The  power  spectrum  (0-100  Hz)  computed  from  the  MVAR  model.  C: 
computed  DTF  values  plotted  against  frequency  (0-100  Hz),  where  the  rows 
correspond  to  the  sending  nodes  and  the  columns  correspond  to  the  receiving 
nodes.  In  each  DTF  plot,  the  red  curves  above  the  black  lines  correspond  to 
frequency  values  with  DTF  values  statistically  significant. 
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Finally,  we  examine  a  case  in  which  the  driving  node  oscillates  in  the  gamma  frequency  range 
and  the  receiving  nodes  are  tuned  to  oscillate  in  the  alpha  band  range.  The  power  spectrum  is 
shown  in  Fig.  6B,  and  in  this  case  we  see  that  node  1  induces  a  large  frequency  component  near 
40  Hz  in  the  receiving  nodes.  This  increase  in  power  is  likely  due  to  the  entrainment  of  many 
neurons  in  the  receiving  node  to  the  driving  node's  frequency.  Examining  the  components  of  the 
DTF  associated  with  these  connections,  we  see  that  they  are  peaked  near  the  50-Hz  region  near 
the  frequency  of  the  driving  oscillating  node.  Analogous  to  the  case  shown  in  Fig.  4  with  nodes 
oscillating  in  the  gamma  frequency  range,  they  are  also  all  substantially  above  the  significance 
level.  There  are  2  prominent,  spurious  components  of  the  DTF  indicating  connections  from 
nodes  4  A  2  and  4  ^  3.  In  addition,  we  see  that  there  are  several  components  of  the  DTF  that  are 
significant  at  high  frequencies  even  though  there  is  no  substantial  power  in  the  spectral  plots. 
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Fig.  6  Gamma  node  1  driving  3  alpha  nodes.  A;  connectivity  structure  of  a 

simulated  network.  B;  power  spectmm  (0-100  Hz)  computed  from  the 
MVAR  model.  C;  computed  DTF  values  plotted  against  frequency 
(0-100  Hz),  where  the  rows  correspond  to  the  sending  nodes  and  the 
columns  correspond  to  the  receiving  nodes.  In  each  DTF  plot,  the  red  curves 
above  the  black  lines  correspond  to  frequency  values  with  DTF  values 
statistically  significant. 
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3.2  Networks  with  Feedback  Loops 

These  2  networks  consisted  of  a  central  driving  node  that  projects  to  3  other  receiving  nodes  with 
feedback  connections  and  a  ring  structure.  Just  like  the  simple  networks,  the  parameters  of  the 
nodes  were  tuned  so  that  the  nodes  would  intrinsically  oscillate  either  in  the  high  alpha  band 
from  10  to  15  Hz  or  in  the  gamma  range  around  40  Hz.  We  examined  the  following  2  cases  with 
a  mismatch  in  frequency  between  the  driving  and  receiving  nodes: 

•  A  gamma  frequency  node  drives  3  alpha  frequency  nodes  with  feedback  and  a  ring 
connection. 

•  An  alpha  frequency  node  drives  3  gamma  frequency  nodes  with  feedback  and  a  ring 
connection. 

In  the  first  simulation,  a  gamma  frequency  node  1  drives  3  alpha  frequency  nodes,  as  shown  in 
Fig.  7,  with  feedback  connections  from  2  of  the  nodes  (3  and  4)  and  an  additional  direct 
connection  between  nodes  2  and  4.  The  power  spectrum  in  Fig.  7B  shows  that  node  1  introduces 
substantial  gamma  frequency  components  in  the  power  spectrum  of  its  3  receiving  nodes.  Thus 
the  components  of  the  DTF  shown  in  Fig.  7C  that  correspond  to  connections  from  node  1  and 
projecting  to  nodes  2-4  are  very  prominent  and  substantially  above  the  level  of  significance  in 
the  frequencies  where  power  exists  in  the  connected  nodes.  However,  just  as  in  Fig.  6,  there  are 
spurious  connections  detected  as  indicated  by  the  significant  DTF  values  at  high  frequencies 
(even  100  Hz)  even  though  there  is  no  substantial  power  in  the  power  spectrum  of  any  of  the 
nodes  at  these  frequencies. 

The  DTF  values  in  Fig.  7C  do  not  robustly  capture  the  feedback  connections  from  3  and  4  to 
node  1.  No  significance  is  found  in  the  connection  3^1,  and  the  connection  4  A  1  is  barely 
detected  by  the  connectivity  measure  and  only  at  the  edge  of  the  higher  frequencies  showing 
power  in  the  spectral  plot.  The  final  connection  2  A  4  shows  significant  DTF  values  at  both 
frequency  ranges  where  power  is  seen  in  the  spectral  plot.  Finally,  there  is  a  spurious  connection 
detected  projecting  from  node  2  to  node  1.  In  this  case,  the  DTF  has  peaks  in  both  the  alpha  and 
the  gamma  range.  These  peaks  correspond  to  the  peaks  of  the  power  spectrum  of  the  sending  and 
receiving  nodes.  Even  though  there  is  a  true  connection  from  node  4  to  node  1,  the  corresponding 
component  of  the  DTF  shows  that  node  4  has  little  influence  on  node  1.  It  is  possible  that  this  is 
an  example  of  the  activity  of  node  2  propagating  around  the  ring  to  node  1 .  However,  this  does 
not  seem  likely  since  the  spurious  peaks  are  more  significant  than  the  4  A  1  component. 
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Fig.  7  Gamma  node  1  driving  3  alpha  nodes  with  feedback  loop.  A:  connectivity  structure 
of  a  simulated  network.  B:  power  spectrum  (0-100  Hz)  computed  from  the  MVAR 
model.  C:  computed  DTF  values  plotted  against  frequency  (0-100  Hz),  where  the 
rows  correspond  to  the  sending  nodes  and  the  columns  correspond  to  the  receiving 
nodes.  In  each  DTF  plot,  the  red  curves  above  the  black  lines  correspond  to 
frequency  values  with  DTF  values  statistically  significant. 
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The  second  simulated  network  in  Fig.  8  has  the  same  connectivity  pattern  as  Fig.  7  but  the 
oscillating  frequencies  of  the  nodes  are  switched:  In  this  case,  node  1  oscillates  in  the  alpha  band 
and  nodes  2-4  oscillate  in  the  gamma  band.  Here  we  see  similar  results  as  those  obtained  in  the 
simple  networks.  Namely,  when  node  1  oscillates  in  the  alpha  band,  it  does  not  entrain  a 
substantial  portion  of  neurons  of  the  higher  oscillating  frequency  nodes.  Thus  there  is 
substantially  less  overlap  in  the  power  spectra  of  the  sending  and  receiving  nodes.  The  resulting 
DTF  components  are  not  significant  (1^3  and  4)  or  only  slightly  above  the  significance  level  at 
the  higher  frequencies  of  oscillation  with  substantial  power  in  node  2(1^2).  However,  the 
feedback  connections  (1  <-  3  and  4)  and  ring  connection  (2  4)  from  the  nodes  oscillating  at  the 

gamma  frequency  are  captured  relatively  well,  although  frequencies  outside  the  peaks  of  power 
in  the  spectral  plot  are  shown  as  significant.  Finally,  we  do  have  some  spurious  connections 
detected  between  (4  <-2)  and  the  one  seen  in  the  first  simulation  in  Fig.  7  between  node  2  to 
node  1.  This  latter  component  of  the  DTF  is  actually  more  significant  than  what  was  observed  in 
Fig.  7. 
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Fig.  8  Alpha  node  1  driving  3  gamma  nodes  with  feedback  loop.  A;  connectivity 

structure  of  a  simulated  network.  B;  power  spectrum  (0-100  Hz)  computed  from 
the  MVAR  model.  C:  computed  DTF  values  plotted  against  frequency 
(0-100  Hz),  where  the  rows  correspond  to  the  sending  nodes  and  the  columns 
correspond  to  the  receiving  nodes.  In  each  DTF  plot,  the  red  curves  above  the 
black  lines  correspond  to  frequency  values  with  DTF  values  statistically 
significant. 
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3.3  Results  Constrained  by  Spectral  Power 

Across  all  6  simulated  networks,  spurious  connections  are  detected;  however,  a  majority  of  these 
spurious  connections  are  detected  at  frequencies  outside  of  the  broadly  tuned  alpha  and  gamma 
oscillations  of  the  network  nodes.  We  propose  that  these  false  positives  may  be  mitigated  by 
examining  the  connectivity  values  only  in  frequencies  that  have  clear  power  peaks  in  the 
spectrum  indicating  presence  of  rhythmic  activity.  For  these  test  networks,  the  maximum  of  the 
power  spectrum  corresponds  to  the  frequency  around  which  most  of  the  cells  in  the  population 
are  synchronized.  Thus,  we  assume  that  the  largest  synchronized  cluster  of  cells  within  a  node 
will  have  the  most  prominent  effect  on  the  receiving  nodes  and  minimal  influence  from  the 
smaller  clusters  oscillating  at  different  frequencies.  Here,  we  constrain  the  analysis  for  each  of 
the  6  simulations  to  the  frequency  range  for  alpha  and  gamma  where  the  power  spectrum  shows 
peaks.  The  range  of  frequencies  for  this  analysis  was  chosen  as  full  width  at  half-maximum  of 
the  powers  density  in  the  node  with  the  highest  values  in  each  frequency  band,  and  we  examine 
how  this  constraint  influences  the  reliability  of  the  DTF  measure. 

Figure  9  shows  the  DTF  results  of  our  6  simulations  for  the  peak  in  the  alpha  range  (10-12  Hz) 
in  light  blue  and  the  gamma  range  (35-45  Hz)  in  dark  blue.  Each  panel  represents  the  results  for 
networks  represented  in  previous  figures  (i.e.,  panel  A  is  data  from  Fig.  3,  panel  B  is  Fig.  4,  etc.). 
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D.  Gamma  node  drives  3  alpha  nodes 
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E.  Gamma  node  drives  3  alpha  nodes  +  feedback 
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F.  Alpha  node  drives  3  gamma  nodes  +  feedback 
_>i  _>2  ~>3  ->4 
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Fig.  9  The  computed  DTF  values  plotted  against  frequency,  where  the  rows  correspond  to  the  sending 
nodes  and  the  columns  correspond  to  the  receiving  nodes.  In  all  of  these  panels,  the  frequency 
bands  used  for  integration  of  significant  DTF  values  were  chosen  as  10-12  Hz  for  alpha  band 
(light  blue)  and  35^5  Hz  for  gamma  band  (dark  blue).  In  all  panels,  node  1  is  the  driving 
node.  The  integrals  in  each  panel  are  labeled  corresponding  to  their  associated  figure;  A  = 

Fig.  3,  B  =  Fig.  4,  C  =  Fig.  5,  D  =  Fig.  6,  E  =  Fig.  7,  and  F  =  Fig.  8. 
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The  simple  networks  are  shown  in  panels  A-D  where  one  driving  node  projects  to  3  other 
receiving  nodes.  The  first  2  networks  in  panels  A  and  B  have  driving  and  receiving  nodes  that 
oscillate  at  the  same  intrinsic  frequency.  As  such,  they  have  strongly  overlapping  spectra,  and 
their  connectivity  pattern  is  captured  relatively  well  by  the  DTP  measure.  While  there  are  some 
spurious  connections  indicated  as  significant  with  Theiler  method,  the  real  connection  from  node 
1  to  all  others  is  clearly  dominant  in  both  integrated  frequency  ranges.  However,  the  same 
network  connectivity  is  shown  in  panels  C  and  D  but  with  a  difference  in  intrinsic  frequency  of 
oscillations  between  driving  and  receiving  nodes,  and  this  configuration  shows  a  different 
connectivity  pattern.  In  panel  C,  the  driving  node  oscillates  at  a  frequency  in  the  alpha  range  and 
the  receiving  nodes  oscillate  in  the  gamma  range.  This  configuration  is  not  detected  well  by 
DTP,  with  mildly  significant  DTP  values  for  only  one  of  the  3  actual  connection  (1  A  3)  and  a 
spurious  connection  (4  ^  2).  In  panel  D,  the  frequency  bands  flip  from  panel  C  with  the  driving 
node  oscillating  in  the  gamma  range  and  receiving  nodes  in  the  alpha  range.  This  configuration  is 
detected  substantially  better  than  in  C,  matching  the  levels  of  significance  seen  with  a  pure 
gamma  network  in  panel  B.  The  real  connection  from  node  1  to  all  others  is  clearly  dominant  in 
both  the  integrated  frequency  ranges  and  the  much  smaller  significance  values  for  the  spurious 
connections  from  (4^1  and  2). 

The  more  complex  networks  with  feedback  connections  and  a  ring  structure  are  shown  in  panels 
E-P.  These  last  2  panels  show  that  significant  DTP  values  do  not  correctly  represent  the 
underlying  connectivity.  In  panel  E,  the  connectivity  pattern  showing  node  1  as  driving  all  others 
is  reproduced  with  higher  DTE  values  in  alpha  band  despite  node  1  dominant  oscillations  in 
gamma  band.  The  feedback  connections  (3  and  4  A  1)  and  connection  2  A  4  are  not  detected,  but 
the  spurious  connection  (2  A  1)  is  captured  with  a  barely  significant  DTE  value  in  the  alpha 
range.  In  panel  E,  when  node  1  is  oscillating  in  alpha  range,  the  connectivity  is  misrepresented  to 
the  largest  degree.  Node  1  is  correctly  shown  to  receive  input  not  only  from  nodes  3  and  4,  but 
also  incorrectly  shown  as  receiving  input  from  node  2.  The  connection  2  A  4  is  correctly 
detected;  however,  the  connections  1  A  2,  3,  4  are  not  detected  at  all. 


4.  Discussion 


In  this  study,  we  simulate  6  networks  using  conductance-based  model  neurons  that  oscillate  in 
frequency  bands  commonly  investigated  in  experimental  EEG  studies,  alpha  and  gamma  (Nunez 
and  Srinivasan  2006).  Eour  of  these  simulated  networks  have  a  simple  feed-forward  structure 
with  one  node  driving  the  other  3,  and  2  have  a  more  complex  connectivity  structure  with 
additional  feedback  connections  and  a  ring  structure.  Across  all  6  simulated  configurations,  we 
manipulate  whether  the  driving  nodes  oscillate  at  alpha  or  gamma.  The  time  series  from  each 
simulated  network  is  separately  fit  to  an  MVAR  model,  and  an  effective  connectivity  measure. 
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DTF,  is  computed  from  the  model  coefficients.  Our  analysis  examined  how  well  this  effective 
connectivity  measure  captured  the  underlying  functional  connectivity  of  the  network  to  assess  its 
utility  for  recovering  network  connections  in  experimental  EEG  data  in  future  research. 

Our  results  demonstrate  that  these  functional  connectivity  measures  can  provide  useful 
information  about  causal  influences  among  network  nodes,  but  they  should  be  applied 
conservatively  and  interpreted  carefully  since  they  detect  spurious  connections  even  in  simplified 
network  configurations.  We  propose  2  guiding  principles  for  application  and  interpretation  of 
these  measures  on  experimental  EEG  data:  1)  analyze  only  those  frequency  bands  with 
substantial  overlap  in  the  power  spectra  of  the  experimental  recordings  from  investigated  nodes 
and  2)  apply  a  conservative  statistical  threshold.  The  measures  performed  best  on  the  simple 
feed-forward  networks,  so  they  are  likely  most  applicable  for  recovering  less  complicated 
network  interactions  and  may  perform  better  when  knowledge  about  anatomical  connections  can 
be  used  to  constrain  their  application  or  interpretation. 

Across  these  6  simulations,  the  measures  correctly  identify  connections  when  the  driving  node 
induces  substantial  activity  in  the  same  frequency  as  the  receiving  node,  resulting  in  significant 
overlap  of  power  spectra  of  both  nodes.  This  effect  was  evidenced  in  both  simple  networks  with 
all  nodes  oscillating  in  the  same  frequency  band  (alpha  in  Eig.  3  and  gamma  in  Eig.  4)  and  when 
the  driving  node  entrained  the  receiving  nodes,  which  was  more  common  for  nodes  oscillating  in 
the  gamma  frequency  (Eigs.  6  and  8).  Eikewise,  in  Eig.  5,  there  was  not  much  overlap  in  the 
alpha  power  spectrum  and  DTE  was  minimally  significant  and  did  not  reliably  capture  the 
connectivity  among  the  nodes.  This  overall  difference  between  alpha  and  gamma  oscillations 
may  heavily  depend  on  the  dynamics  of  the  networks  in  question.  In  our  model,  afferent 
connections  from  a  network  oscillating  at  a  high  frequency  to  a  network  oscillating  at  a  lower 
frequency  are  detected  at  a  much  weaker  coupling  strength  than  in  the  opposite  scenario.  Thus, 
our  model  networks  have  a  preference  for  higher  frequency  inputs,  and  this  highlights  a  critical 
element  of  our  results:  the  importance  of  power  overlap  in  the  spectrum. 

Our  results  show  that  entrainment  between  the  nodes  is  not  limited  to  their  dominant  oscillating 
frequency.  The  measure  also  detected  connectivity  in  other  frequency  ranges  since  the  driving 
node  may  synchronize  over  the  total  activity  in  the  receiving  node  so  the  effective  connectivity  is 
visible  across  the  whole  frequency  range.  Consequently,  we  propose  that  the  analysis  should  be 
constrained  only  to  those  frequencies  with  overlapping  power  in  the  spectrum.  This  enhanced  the 
level  of  significance  for  the  connectivity  values,  as  demonstrated  in  Eig.  9. 

Of  course  limiting  the  analysis  to  overlapping  power  in  the  spectrum  may  miss  some  real 
connections.  In  an  experimental  data  set  there  may  be  synchronous  activity  in  a  frequency  band 
not  corresponding  to  the  maximum  of  the  power  spectrum  (Brovelli  et  al.  2004).  In  this  case,  the 
ability  of  one  population  of  cells  to  entrain  another  at  a  certain  frequency  is  dependent  both  on 
anatomical  connectivity  as  well  as  the  biological  properties  of  the  cells  involved  in  the 
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oscillation.  In  such  cases,  one  would  have  to  employ  a  variety  of  information  to  attempt  to  infer 
the  effective  connections,  including  any  details  about  the  true  anatomical  connectivity  between 
cell  populations  and  brain  regions. 

Although  implementing  a  constraint  based  on  the  power  spectrum  improves  the  performance  of 
the  connectivity  measures,  using  a  statistical  method  with  significance  level  0.05  still  shows 
some  spurious  connections.  In  more  complicated  networks,  when  subnetworks  at  each  node  have 
different  levels  of  synchrony  or  oscillate  at  different  frequencies,  the  measures  may  function  less 
reliably.  This  was  evident  in  Figs.  7  and  8,  which  indicated  an  influence  from  node  2  to  node  1 
where  no  connection  exists.  To  maximize  the  correct  number  of  detected  connections,  one  must 
not  only  examine  the  measures  in  a  narrow  band  around  the  peak  of  the  power  spectrum,  but  also 
set  a  relatively  conservative  significance  level.  Future  work  should  examine  if  alternative 
statistical  approaches  or  surrogate  data  approaches  can  improve  identification  of  spurious 
connections. 

When  interpreting  connectivity  analyses  of  experimental  data,  these  connectivity  measures  are 
likely  best  used  to  confirm  a  hypothesis  about  effective  connectivity  based  on  some  other 
anatomical  and  physiological  data.  These  measures  are  not  well  suited  to  be  used  in  isolation  to 
try  and  determine  unknown  anatomical  connections,  but  they  can  likely  reveal  dynamic  casual 
interactions  between  brain  regions  when  there  is  already  some  knowledge  to  formulate 
hypotheses  about  the  underlying  network  connectivity.  Previous  research  has  shown  that 
bivariate  Granger  Causality  fails  to  correctly  detect  functional  connections  in  a  simple  Pyloric 
central  pattern  generator  (Kispersky  et  al.  2011).  They  computed  the  Granger  Causality  measure 
on  a  time  series  obtained  from  the  stomato-gastric  ganglion  of  the  crab,  which  consists  of  3 
oscillating  subpopulations  of  neurons.  They  found  that  the  computed  Granger  Causality  inferred 
casual  influences  where  no  anatomical  connections  exist  (Kispersky  et  al.  2011).  Thus,  our 
results  are  complementary  and  suggest  that  methods  based  on  autoregressive  modeling  can 
produce  spurious  results  when  feedback  connections  are  introduced  into  a  network  or  when  the 
network  is  more  complicated  than  a  simple  network  with  a  single  source  projecting  to  other 
nodes. 

Furthermore,  these  linear  methods  are  prone  to  error  simply  from  the  fact  that  we  are  attempting 
to  fit  a  time  series  generated  from  a  nonlinear  process  to  a  linear  model.  Thus,  this  may  be 
another  reason  that  we  detect  spurious  connections.  In  particular,  DTF  and  similar  methods 
reliant  on  computations  in  the  frequency  domain  are  unable  to  detect  cross-frequency 
interactions,  but  in  nonlinear  systems  like  the  brain,  there  might  be  interactions  between  systems 
not  inducing  the  same  frequency  activity  in  the  receiving  node.  Future  research  should  examine 
nonlinear  alternatives  for  recovering  functional  connectivity  patterns. 
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5.  Conclusion 


Our  simulations  suggest  the  directed  connectivity  measures  should  be  interpreted  carefully  and 
applied  conservatively  to  experimental  EEG  datasets.  We  suggest  2  guiding  principles  for 
experimental  analyses.  Eirst,  DTE  should  only  be  analyzed  in  a  frequency  band  centered  at  the 
frequency  ranges  with  substantial  overlap  in  the  power  spectra  of  the  signals  from  different  brain 
regions.  Second,  a  conservative  estimate  of  the  significance  should  be  used  to  avoid  detecting 
spurious  connections. 

Although  these  guiding  principles  will  not  completely  eliminate  spurious  results,  they  may 
greatly  reduce  the  number  of  false  positives.  These  measures  may  be  particularly  helpful  when 
testing  hypotheses  about  connectivity  among  known  anatomical  structures,  but  they  are  unlikely 
to  be  informative  when  trying  to  infer  underlying,  unknown  connectivity.  Additional  measures 
such  as  phase-synchronization  measures  or  measures  based  on  autoregressive  models  with  time- 
dependent  coefficients  may  be  further  used  as  a  check  on  the  validity  of  the  results  obtained  with 
measures  similar  to  DTE.  The  inclusion  of  such  measures  may  provide  yet  another  way  to  bolster 
our  criterion  for  validating  whether  or  not  2  brain  regions  are  interacting. 
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The  following  are  the  values  of  the  parameters  used  in  the  neuron  model  based  on  the  model 
developed  by  Kudela  et  al.'’^  The  individual  channel  currents  are  described  by  the  following 
equations: 

Inu  =  gNamlVdl  -  W)(V  -  Vr,a), 

Ica^gCaX^  ^(V-Vca), 

lK  =  gKW\V-V„),  (1) 

lA^gAAc.(V)B(V-V„), 
and 

iL=gL(y-v,). 

The  recovery  and  activation  variables  for  the  previous  equations  are  described  by 

dF  Foo(V)-F 

—  =  ,  (2) 

ut  Tp 

where  F  represents  W,  X,  A,  m,  and  B.  The  steady  state  functions  for  these  variables  are  given  by 
an  expression  of  the  form 

(3) 

l+exp\^-2af'{v-v[/2)) 

In  all  equations,  V  is  the  neural  membrane  potential. 

The  synaptic  current  Isyn  is  modeled  as 

4yn(t)  =  Wjgj{t  -  Tj)  {Vrest  “  E^y^),  (4) 

where 

gjit)  =  gsyn  -  t)/Td)  -  exp((tfc  -  O/To)),  (5) 

where  k  denotes  summation  over  past  action  potentials  on  a  given  synapse  and  j  denotes 
summation  over  all  synaptic  inputs  of  the  neuron. 

The  values  of  conductances  and  reversal  potentials  used  in  the  model  neuron  are  as  in 
Kudela  et  al.^ 


*  Kudela  P,  Franaszczuk  PJ,  Bergey  GK.  A  simple  computer  model  of  excitable  synaptically  connected  neurons.  Biological 
Cybernetics.  1997;77(l):71-77. 

2 

Kudela  P,  Franaszczuk  PJ,  Bergey  GK.  Changing  excitation  and  inhibition  in  simulated  neural  networks:  effects  on  induced 
bursting  behavior.  Biological  Cybernetics.  2003;88(4):276-285.  doi:10.1007/s00422-002-0381-7. 
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In  Table  A-1  are  listed  the  values  that  specify  the  number  and  types  of  connections  between  the 
excitatory  and  inhibitory  neural  subnetworks.  Node  connections  are  listed  in  Table  A-2. 

Table  A-1  Parameter  values  for  normal  alpha  activity  in  each  subnetwork 


Connection  Type 

Domain  Size 

Connections  Strength 

STD  Delay  STD 

Excitatory  to  excitatory 

4x4 

30  5  3  100  50 

Excitatory  to  inhibitory 

4x4 

30  20  10  100  50 

Inhibitory  to  excitatory 

4x4 

30-100  50  100  50 

Inhibitory  to  inhibitory 

4x4 

5-100  75  100  50 

Note:  STD  is  a  standard  deviation  of  Gaussian  distribution  of  the  parameters. 


Table  A-2  Connections  between  nodes 


Connection  Type 

Domain  Size 

Connections  Strength 

STD  Delay  STD 

Excitatory  to  excitatory 

1  X  1 

30  100  no  150  50 

Excitatory  to  inhibitory 

1  X  1 

10  100  no  150  50 

Inhibitory  to  excitatory 

1  X  1 

10-60  no  150  50 

Inhibitory  to  inhibitory 

1  X  1 

20-60  no  150  50 
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(1) 


The  formulas  for  MVAR  model  order  selection  criteria  used  in  this  study  are  the  following: 

AIC(p)  =  log{det(V))  +  2^, 


fi/C(p)  =  log{det(y))  + 
HQICip)  =  log(dem)  + 


and 


log{FPE(p))  =  (op(det(K(p)))  + 
AICc(p)  =  AIC(p)  +  2 

N-K(p)-1 


BICc{p)  =  fi/C(p)  + 


log(w-A:(p)-i)ft:(p)(A:(p)+i) 

(N-K(p)-l) 


HQICcip)  =  HQICcip)  + 


\og{log(,N-K{p)-l)')K{p){K{p)+l) 

iN-K(p)~l) 


(2) 


(3) 

(4) 

(5) 

(6) 


In  these  equations,  K(p)  =  pM^+M(M+l)/2  is  a  number  of  coefficients  estimated  for  the  model 
order  p  for  M  channels,  and  N  is  the  number  of  sample  points  used  in  estimation.  All  of  the 
model  order  selection  criteria  are  derived  under  the  assumption  that  the  data  in  the  time  series 
x(n)  is  Gaussian  distributed.^ 


Stoica  P,  Selen  Y.  Model-order  selection:  a  review  of  information  criterion  rules.  IEEE  Signal  Processing  Magazine. 
2004;21(4):36^7.  doi:10.1109/MSP.2004.1311138. 
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